About Forum Issues Tutorials Documentation </span>

Example 3: Simulating a Holliday Junction PDB assembly

This notebook takes a crystal structure from the PDB and prepares it for simulation.

Author: Aaron Virshup, Autodesk Research
Created on: July 1, 2016
Tags: DNA, holliday junction, assembly, PDB, MD



In [ ]:

    
%matplotlib inline
from matplotlib.pyplot import *

import moldesign as mdt
from moldesign import units as u

A. View the crystal structure
B. Build the biomolecular assembly
C. Isolate the DNA
D. Prep for simulation
E. Dynamics - equilibration
F. Dynamics - production

A. View the crystal structure

We start by downloading the 1KBU crystal structure.

It will generate several warnings. Especially note that it contains biomolecular "assembly" information. This means that the file from PDB doesn't contain the complete structure, but we can generate the missing parts using symmetry operations.



In [ ]:

    
xtal = mdt.from_pdb('1kbu')
xtal.draw()

B. Build the biomolecular assembly

As you can read in the warning, 1KBU only has one biomolecular assembly, conveniently named '1'. This cell builds and views it:



In [ ]:

    
assembly = mdt.build_assembly(xtal,'1')
assembly.draw()

By evaulating the "assembly" object (it's a normal instance of the moldesign.Molecule class), we can get some information about it's content:



In [ ]:

    
assembly

Because we're only interested in DNA, we'll create a new molecule using only the DNA residues, and then assign a forcefield to it.

C. Isolate the DNA

This example will focus only on the DNA components of this structure, so we'll isolate the DNA atoms and create a new molecule from them.

We could do this with a list comprehension, e.g. mdt.Molecule([atom for atom in assembly.atoms if atom.residue.type == 'dna'])

Here, however we'll use a shortcut for this - the molecule.get_atoms method, which allows you to run queries on the atoms:



In [ ]:

    
query = {'residue.type':'dna'}
dna_atoms = assembly.get_atoms(**query)
dna_only = mdt.Molecule(dna_atoms)
dna_only.draw3d(display=True)
dna_only

D. Prep for simulation

Next, we'll assign a forcefield and energy model, then minimize the structure.



In [ ]:

    
dna = mdt.assign_forcefield(dna_only)



In [ ]:

    
dna.set_energy_model(mdt.models.OpenMMPotential, implicit_solvent='obc')
dna.configure_methods()



In [ ]:

    
minimization = dna.minimize()



In [ ]:

    
minimization.draw()

E. Dynamics - equilibration

The structure is ready. We'll associate an integrator with the molecule, then do a 2 step equilibration - first freezing the peptide backbone and running 300K dynamics, then unfreezing and continuing dyanmics.



In [ ]:

    
# Freeze the backbone:
for residue in dna.residues:
    for atom in residue.backbone:
        dna.constrain_atom(atom)



In [ ]:

    
dna.set_integrator(mdt.integrators.OpenMMLangevin,
                   timestep=2.0*u.fs,
                   frame_interval=1.0*u.ps,
                   remove_rotation=True)
dna.integrator.configure()

And now we run it. This is may take a while, depending on your hardware.



In [ ]:

    
equil1 = dna.run(20.0*u.ps)



In [ ]:

    
equil1.draw()

Next, we'll remove the constraints and do full dynamics:



In [ ]:

    
dna.clear_constraints()
equil2 = dna.run(20.0*u.ps)



In [ ]:

    
equil = equil1 + equil2
equil.draw()



In [ ]:

    
plot(equil2.time, equil2.rmsd())
xlabel('time / fs'); ylabel(u'rmsd / Å'); grid()

NOTE: THIS IS NOT A SUFFICIENT EQUILIBRATION FOR PRODUCTION MOLECULAR DYNAMICS!

In practice, before going to "production", we would at least want to run dynamics until the RMSD and thermodynamic observabled have converged. A variety of equilibration protocols are used in practice, including slow heating, reduced coupling, multiple constraints, etc.

F. Dynamics - production

Assuming that we're satisfied with our system's equilibration, we now gather data for "production". This will take a while.



In [ ]:

    
trajectory = dna.run(40.0*u.ps)



In [ ]:

    
trajectory.draw()

G. Save your results

Any MDT object can be saved to disk. We recommend saving objects with the "Pickle" format to make sure that all the data is preserved.

This cell saves the final trajectory to disk as a compressed pickle file:



In [ ]:

    
trajectory.write('holliday_traj.P.gz')

To load the saved object, use:



In [ ]:

    
traj = mdt.read('holliday_traj.P.gz')



In [ ]:

    
traj.draw()